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Basic principles of hp Virtual Elements on quasiuniform 

meshes 

L. Beirao da Veiga * * * § A. Chernov } L. Mascotto f A. Russo ® 


Abstract 

In the present paper we initiate the study of hp Virtual Elements. We focus on the case 
with uniform polynomial degree across the mesh and derive theoretical convergence estimates 
that are explicit both in the mesh size h and in the polynomial degree p in the case of finite 
Sobolev regularity. Exponential convergence is proved in the case of analytic solutions. The 
theoretical convergence results are validated in numerical experiments. Finally, an initial study 
on the possible choice of local basis functions is included. 


1 Introduction 

The Virtual Element Method (VEM) is a very recent generalization of the Finite Element Method, 
introduced in j5], that responds to the increasing interest in using general polyhedral and polygo¬ 
nal meshes, also including non convex elements and hanging nodes. The main idea of VEM is to 
use richer local approximation spaces that include (but are typically not restricted to) polynomial 
functions and, most importantly, avoid the explicit integration of the associated shape functions. 
Indeed, the operators and matrices appearing in the problem are evaluated by introducing an in¬ 
novative construction that only requires an implicit knowledge of the local shape functions. By 
following such developments, the VEM acquires very interesting properties and advantages with 
respect to more standard Galerkin methods, yet still keeping the same implementation complex¬ 
ity. For instance, in addition to allowing for polygonal and polyhedral meshes, it can handle 
approximation spaces of arbitrary C k global regularity on unstructured meshes. 

Although the Virtual Element Method has been applied to a large range of problems (a non 
exhaustive list being ffil4ll MT^irMT7l[2Rl22ll25l 1. all the present works on VEM are focused, both 
theoretically and numerically, on the h —behaviour of the method. In other words, the convergence 
properties of the schemes are investigated assuming that the polynomial degree p is fixed and only 
the mesh is refined. On the other hand, looking at the Finite Element literature, a very successful 
approach in applications is to allow for a variable value of p and to focus not only in the accuracy 
that can be obtained by reducing the mesh size h , but also by increasing p. As in the FEM 
literature, we here refer to such approach as hp analysis; we mention for instance 
as very short list of papers and books among the very large literature of hp FEM. 

The aim of the present paper is to initiate the study of hp Virtual Element Methods. The first 
motivation of such study is to show that the powerful hp methodology can be adopted also in the 
framework of Virtual Elements. The second, but not secondary, motivation is that we believe that 
combining the huge mesh flexibility of VEM with the advantages of a full (possibly adaptive) hp 
method can yield a very efficient and competitive methodology. 

The present contribution focuses on the initial foundations of such ambitious plan, mainly in 
terms of convergence estimates. We restrict our attention to a two dimensional scalar elliptic 
model problem (as in m) and assume a polynomial degree p that is the same on all elements of the 
(quasi uniform) mesh. First of all, we prove fundamental convergence results (and the associated 
interpolation estimates) that is explicit in both h and p. As a second result, we show that also 
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Virtual Elements can attain exponential convergence when the target solution is analytic on a 
suitable (small) extension of the domain. We then explore numerically the behaviour of Virtual 
Elements in terms of p, both in the case of solutions with finite Sobolev regularity and for analytic 
solutions. In the Appendix, we start to explore another interesting issue of hp elements, that is 
the choice of the basis and the condition number of the associated stiffness matrix. Note that, 
since in this work we focus on scalar problems in a planar two-dimensional domain, direct solvers 
can generally be used and the condition number issue is not of primary importance. Indeed, it 
is mainly the stability of the solver that determines the best attainable accuracy, as we show in 
the numerical tests. Nevertheless, in order to answer to some natural questions (such as: how do 
Legendre bases cope on general polygons?) we decided to include an initial study related to the 
choice of the basis. 

The paper is organized as follows. After presenting the continuous model problem, in Section 
[3] we make a brief review of the Virtual Element Method. Afterwards, in Section 0 ] we present the 
theoretical error estimates and in Section[G]we develop the associated numerical tests. Finally, the 
Appendix follows. 


2 The model problem 

Let O be a simply connected polygonal domain and let T be its boundary. Let H 1 (oj ), with l £ No 
and cj open measurable set, denote the usual Sobolev space with square integrable weak derivatives 
of order /; let || ■ and | ■ denote the associated norm and seminorm, respectively (see [T]). 
Let / £ L 2 (fl). We consider the two dimensional Poisson problem with homogeneous Dirichlet 
boundary conditions: 

— A u = f in Cl, it = 0 on r. (1) 

We set V := Hq(CI) and we consider the weak formulation of problem Jl]) 

find u£V such that a(u,v) = (f,v) o,n, Vu £ V, (2) 

where (-,-)o,n is the L 2 scalar product on and a(-,-) := (V-,V-)o,n- 

It is well known that problem @ is well-posed (see for instance m) since the bilinear form a 
is continuous and coercive (i.e. a(v,v) > a||u||f q, where a > 0) thanks to the Poincare inequality. 
Throughout this paper, C denotes a positive constant whose dependence on certain parameters 
will be made explicit where necessary. 


3 Virtual Elements for the Poisson problem 

In the present section, we introduce a Virtual Element Method for the Poisson problem @ based on 
polygonal meshes. Let (7 ~h}h be a sequence of decompositions of Cl into non-overlapping polygonal 
elements K of diameter hx = diam(A) := sup{|x — y| : x, y £ K}. The characteristic mesh size 
is denoted by h := ma x{hx ■ K £%,}■ Let Vh and £h be the sets of all vertices and edges in the 
mesh Th respectively. Moreover, we denote by := Vh D dCl be the set of all boundary vertices 
and by the set of edges e of an element I\ £Th- 

Henceforth, we assume that there exist two positive real numbers 7 and 7 such that the sequence 
of decompositions satisfies the following: 

(DO) the decomposition Th is made of a finite number of simple polygons of diameter hx , 

(Dl) for all K £ Th, K is star-shaped with respect to a ball of radius > hx"f, 

(D2) for all K £ Th, the distance between any two vertices of K is > hxl- 

To every edge e £ £h we associate a tangential vector r e and a normal unit vector n e obtained by 
a counter-clockwise rotation of T e . 

We split the bilinear form a as a sum of local contributions 

a(u,v ) := a K (u,v), Vu, v £ V, 

KeT h 
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with a K (u,v) := (Vit, Vv)q,k■ 

It was shown in [S] that it is possible to build: 

• V h (K), a finite dimensional subspace of Hq(SY)\k', 

• a symmetric local bilinear form aj£ : 14(A) x 14 (A) —> R; 

• Vh a finite dimensional subspace of Ag(fl) such that Vh\x = V h (K)- 

• a symmetric bilinear form ah ■ Vh x Vh —t R, of the form ah(uh,Vh) = a h ( u hi v h), 

\/uh,Vh £ Vh where a^ : Vh\x x Vh\K —> R are local symmetric bilinear forms; 

• an element fh € V)( and a duality pairing (•; -)h 
in such a way that the resulting discrete problem 

find Uh G Vh such that a h (u h , v h ) = {f h \ v h )h, Vv h £ V h (3) 

has a unique solution Uh € 14 which is close to the solution u of the original problem d2J) ■ More 
precisely, when u £ H k (Q) the error in the energy norm admits the upper bound 

if u £ H k (Q), k > 1 , \u-u h \i } n < C7i fc-1 |u| fci n (4) 

where the constant C = C(p) depends implicitly on the (fixed) polynomial degree p but not on 
the characteristic mesh size h. We now briefly review the local Virtual Spaces introduced in [jjj. 
Let K £ Th and let p £ N, p > 1. Let P p (e) and P p _2 (K) be respectively the set of polynomials of 
degree p over the edge e and of degree p— 2 over the polygon K, with the convention P_i(If) = {0}. 
With the space of contunuous piecewise polynomials over the boundary of each element K 

M p (dK) := {v h £ C°(dK) \ v h \ e £ P p (e), Ve £ } , 

we define the local Virtual Element spaces 

V h (K) := {v h £ H\K) | A v h £ P P - 2 (K), v h \ dK £ B p (5I^)} . (5) 

Observe that P P (K) C Vh(K) for any K £ Th- For any fixed function v £ Vh{K) we identify the 
following set of local degrees of freedom: 

• the values of v at vertices of K , 

• the values of u at (p — 1 ) internal nodes of each edge e £ £ p (for instance at the internal 
Gaufi-Lobatto nodes, as done in 0 and mi), 

• the internal moments f K q a Vhdx, where {q a : 0 < |a| < p(p — l)/ 2 } is abasis for P p _ 2 (A'). 

For instance, [9] and m employed a basis of shifted and scaled monomials: let xk and hx 
be the barycenter and the diameter of K respectively, then g a (x) := for any 

a = (ai, a 2 ) £ No suc h that |a| := a\ + a 2 < p — 2 . 

The global Virtual Space is obtained by the continuous matching of the local spaces over the 
element boundaries 

Vh := {vh £ C°(fi) : Vh\x £ Vh(K),Vh\dn = 0} C Hq(£1) 

with the natural definition of the global degrees of freedom from the local ones. It was shown in [5] 
that, given K £ Th, the bilinear forms a must satisfy the two following assumptions. 

(Al) p-consistency: VFT £ Th it holds that 

a K (q, v h ) = a%(q,v h ) 7 Vv h £ V h (K), \/q £ P P (K); (6) 

(A 2 ) stability: MK £ Th there exist two constants 0 < a* < a* < oo such that 

a*a K (v h ,v h ) < a,h(v h ,v h ) < a*a K {v h ,v h ), Vv h £V h (K). (7) 
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Let ip be a sufficiently regular function, e.g. ip G H 1 {K). We introduce the local averaging operator: 


Wk\ fan Vip)dx, if p = 1 
jk\ Ik v{x)dx, if P > 1 ’ 


( 8 ) 


Having this, we introduce the projection operator nj : Vh(K) P P (K) as follows: for any 

Vh G Vh(K) its projection n Jvh G P P {K) is the unique polynomial satisfying, 

a K (nJv h - v h ,q) = 0 VgGPp(iL), 

- I—V (^) 

nj Vh — Vh = 0 where the averaging operator is defined in © . 

Then, a candidate bilinear form ah satisfying (Al) and (A2) can be sought in the form: 

ah(u h ,v h ) = a K (nj u h ,Ujv h ) + S K (u h - Uju h ,v h - TLjv h ), Vu h ,v h G V h (K), 

where S K is a positive definite bilinear form satisfying 

c 0 a K (v h ,v h ) < S K (v h ,v h ) < C\a K {vh, v h ), Vv h G V h (K), such that Tljv h = 0, (10) 

for some positive constants cq and ci independent on h , p and K. 

The global discrete bilinear form reads 


a h (u h ,v h ) := ^ a h( u h,Vh), Vith, v h € V h - 

K&T h 


Finally, we recall from [jj a possible choice for the loading term. Let Pp'I\ and Pq' K be the 
L 2 -projector on polynomials of degree p — 2 and 0 respectively over the polygon K and let the 
averaging operator be defined in (O. Then, we may define 


{fhi Vh.)h 


J KGT, 
KKGT, 



v h dx, 


v h dx, 


Vv h G V h if p > 2 , 
\/v h G V h if p = 1 • 


( 11 ) 


Under the above choices for Vh, ah and fh, [1] guarantees well-posedness and h -convergence ©• 

Remark 1. As shown in [9], the projection operator nj in d9]) is computable using the degree of 
freedom values, without the need of any further information on the virtual shape functions. We 
finally note that the definition in © is not the only possible one; other (computable) choices could 
be used instead. 


4 Approximation results 

In this section, we give a convergence result for the error of the Virtual Element Method measured 
in the energy norm in terms of both h and p. 


4.1 Auxiliary approximation results 

Let u and Uh be the solutions of © and © respectively, and denote by S^~ l (Th) the space of 
the piecewise discontinuous polynomials of degree p over the decomposition Th- Given u G iL 1 (/\), 
\/K G Th , we define the broken ULseminorm as 

Mm, a = E (Mi,*)*- (12) 

K&T h 

Let Th be the smallest constant satisfying 


(f,Vh)o,n — {fh , Vh)h < Th\vh\i,n, Vvh G Vh- 


(13) 
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Then the following best approximation estimate holds (see Theorem 3.1 in 0) 

\u — Uh\i.ti < C [ inf \u — UttIi h.fi + inf \u - + T h , (14) 

yuvesl’-'in) J 

where C is a constant depending only on ct* and a* from assumption (A2). In what follows, we 
shall derive estimates for the three terms in m that are explicit in both h and p. 


4.1.1 Polynomial approximation term 


We start by bounding the term |u — u^h,i,o.- In order to derive the bound, we need to prove 
a generalized-polygonal version of a classic result, namely Lemma 4.1 in [6\ In this lemma, it 
was shown the existence of a sequence of polynomials which approximate H k functions over the 
triangular and square reference elements. We extend this result for generic polygons having the 
unit diameter. Thus, we are ready to show the following 

Lemma 4.1. Let K C M 2 be a polygon with diam(K)=1. Moreover, assume that K is star-shaped 
with respect to a ball of radius > 7 and the distance between any two vertices of K is >7, 7 and 
7 being the constants introduced in assumptions (Dl) and (D2) of Section^ Then, there exists a 
family of projection operators {11^ }, p = 1, 2,... with 11^ : H k+l (K) —» P P (K) such that, for 
any 0<£<k + Y,uG H k+1 (K), k GN, it holds 




e,K 


< Cp 


~(k+\-l 


\k+l,K 


(15) 


with C a constant independent on u and p. 

Proof. We assume without loss of generality that xg, the barycenter of K , coincides with the 
origin 0. For a given r > 0, we define 


R(r) := {{x,y) G R 2 | \x\ < r, \y\ < r} . 


(16) 


Thanks to the fact that diam(A')=l and x^ = 0, we have R( 1) D K. 

Let ro > 1. Then, it obviously holds K C R(r 0 ). We note that dK is Lipschitz; consequently, 
using EH, there exists E : H k+1 (K) —> H k+1 (R(2ro)) extension operator such that Eu = 0 on 
R(2rg) \ A(|ro) and \\Eu\\k + i R( 2 r 0 ) < C||^||fc + i A careful inspection of Theorem 5 in Chapter 
VI of [2Zj shows that the constant C depends only on k and on the “worst angle” value 


9t> = min minld, 27r — 0|, 
K O^k 


where A^ denotes the set of the (amplitude of) internal angles of K. In particular, the constant C 
may explode when 9^ —> 0. It is easy to check that, under the regularity hypotheses on K , the angle 
parameter 9^ is bounded from below by a constant depending only on 7 , 7 . Thus the constant C 
can be bounded in terms of k and 7 , 7 . Therefore, it holds \\Eu\\ k+1 R ( 2 r 0 ) < C(/C; 7 ; 7 )|| : «||fc_|_i 
The remaining part of the proof, that is based on the approximation of the extended function Eu, 
follows exactly the same steps as in )b], Lemma 4.1, and is therefore not shown. □ 

Using this result, we are able to give a generalized-polygonal version of Lemma 4.5 of [ 6 ], which will 
play the role of local hp estimate result on \u — m w |i k , where AT is a polygon of the decomposition 

%. 

Lemma 4.2. Let K G 7h satisfying assumptions (Dl) and (D2) and u G H k+1 {K). Then there 
exists a sequence of projection operators {11^^}, p = 1,2,..., with : H k+1 (K) —> Pp(A') such 
that for any 0 < i < k + 1, k G N 

\ u ~ n k, p u U,k < G p % +1 _ e |M|fc+i,JO ( 17 ) 

where p, = min(p, k) and C is independent on u, h and p. 
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Proof. We consider the mapping F(x) = ^(x — xk). Let K = F(K), where Hk denotes the 
barycenter of K. Obviously diam(it')=l and the barycenter of K is in the origin, xg = 0. Let 
Ilg p u be the sequence of approximating polynomials of degree p, introduced in Lemma 14.11 We 
set II k p u be the push forward of the above sequence with respect to the transformation F, i.e. 
II ^ p u = (Ilg (u)) o F, where ip = ip o F —1 for a sufficiently regular function p. Then, it is easy 
to check, by a simple change of variables argument, that 

I u ~ ^K,p u \^,K ^ Ch K I U — Lt K,p U \e,K’ 


where C is a constant independent on K (hence on K ), h, u and p; besides, C is independent also 
on £, 7 and 7 , thanks to the fact that F is the composition of a translation with a dilatation. 

We apply Lemma FTTI and we obtain, by adding and subtracting any q € F p (K), 


|u-n 


K,p L 


e,K < Ch K \\(u-q)-U^ p (u-q)\\ e S; < C{k , 7 , 7 ) 


h 1 ~ i 

n K 

n k-\-l—£ 


‘-511 


fe+l,K’ 


\/q€F p (K), 


(18) 

where C in the right hand side of (flSl) is a constant depending on k, 7 , 7 , since clearly K = F(K) 
still satisfies the shape regularity assumptions with the same constants 7 , 7 . Using the classical 
Scott-Dupont theory (see e.g. [15] or [2U]) and a scaling argument, bound (1T51) yields 


l U “ n K,P U kK < 



<c^- 


— p k+1 - ( - 


\k+i,K, /r = min (p,fc), (19) 


where C is independent on u,p and h. 


□ 


Remark 2. We note that if k < p then the classical Bramble-Hilbert Lemma allows to take the 
seminorm in the right hand side of m, yielding 


i — ^ 

\u- ^ h K,pAt,K < C{k, 7 , 7 ) % +1 _ e \u\k+l,K- 

We are now able to give a global estimate on \u — u n \h,i,n in dH, where € S%’ (7ft), 

3%’- (7ft) being defined at the beginning of Section IXTl In fact, by choosing u v \k = P u f° r 
K £ 7ft and recalling the shape regularity properties (D1)-(D2), we obtain: 


\u-u v \ h p,n < Ci-r||«|U+i,n, p = min(p, fc), 

pK 

h k 

\u ~ W 7 r|ft,i,n < C 2 -rMfc+i,n, for p > k 

pK 

where C\ and C 2 are two constants independent on u , p and h. 


( 20 ) 


4.1.2 Virtual Interpolation term 


We turn now to the term |u — in (TUH) . Preliminarily, we observe that (Dl) and (D2), defined 

in Section [3J imply that there exists 7ft, an auxiliary conformal triangular mesh that refines 7ft, 
obtained by connecting, for all K £ 7ft, the Nk vertices to the center of the ball that realizes 
assumption (Dl) for K. Moreover, it is easy to check that each triangle T £ 7ft is uniformly shape 
regular. 

Let 5ft’ (7ft) be the set of continuous piecewise polynomials of degree p over the auxiliary 
triangular decomposition introduced above. It is well known (see Theorem 4.6 in [Bj) that there 
exists p p € S^’ 0 (Th) such that for any u € F[ k+1 (tt), k € K 


1 

||u- <K||i,n < C'l-rlHIfc+i.Q, k>~, 

F p K 2 

h k 1 

\u - Pp\i,n < C 2 —j-\u\ k +i,ci, k > —, p > k 

y p K z 


( 21 ) 
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where C\ and C 2 are two constants independent on u , p and h and where p = min(p, k). 

Now, we use in m to construct an approximant ui £ 14 of u. For this purpose, we modify 
a particular technique introduced in [25] . 

Lemma 4.3. Under (Dl) and (D2), for all u £ H k+1 (Q), k £ N, there exists ui £ 14 such that 


|U - < C—r 

pK 


|u]|fc+i,n, 




( 22 ) 


where C is independent on h, p and u. 

Proof. Let u w be the function defined in (121)1) . Let ip k be the function described in (l2ll> . For each 
K £ 7h, we define ui\k the solution of the following problem: 


{ —A ui = -Au t in K 
ui = ip k on dK 

It is easy to check that ui\k £ 14,| k ■ Moreover, since ui £ 74 1 (fl), it holds that ui £ 14. 
Using (l23l) . we can write 

f —A (uj — u n ) = 0 in K 
( ui — u % = <Pp — Uk on dK 

Therefore, since (ui — u is harmonic it holds 

|Ui - UTr\i,K = inf {\z\i,k, 2 £ H 1 ^) \ z = (fip - Uir on dK) < \ip k - u^k ■ 
Finally by (1211) we obtain 

| U - Ui\i^ K < | U- U^i'K + I^tt - Ui\i,K < \u - U n \ lt K + |«7r - Pp\l,K 
< 2| U - U n \ ltK + \u~ Pp\l,K- 

The proof is completed by summing on all the elements in (1251) and using (1201) , (1211) . 


(23) 


(24) 


(25) 

□ 


Remark 3. We point out that if k < p and under the hypothesis of Lemma 1Q1 the following holds: 

h k 

| u- ui |i,n < C'-r-Mfc+i.fi, 

pK 

C depending only on £:, 7 and 7 (introduced in (Dl) and in (D2)). 


4.1.3 Loading approximation term 

It remains to estimate the term Th in (fill) . We have the following result. 

Lemma 4.4. Under assumptions (Dl) and (D2), let the loading term f £ H k+1 (K) for all 
K £ Th, k £ N. Then it holds 

Ph < C-|— ( II/H|k ) ’ M = min(p,fc + 2). (26) 

P \KeTh ’ ) 

where C is a constant independent on h, p and u. 

Proof. Since the case p = 1 has been already analysed in [5], we only consider the case p > 2. Let 
Vh £ 14. Let Pp’ K 2 be the L 2 -projector on polynomials of degree p — 2 over the polygon K , for all 
K £ Th- We get by (HU) 

{f, v h )o,n ~ ( fh, v h )h = ^2 (f ~ P pl2fi Vh)a,K = (/ - Ppl^f, v h - P°^v h ) 0 ,K 

KeTh KeT h 

< 11/~ Pp’-2f\\o,K\\ V h - Pp-2 v h\\o,K 

KGT h 

< Wf ~ fp-2\K\\o,K\\vh - Vp_ 2 \ K \\o,K, 

K&T h 
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where fp_ 2 \K and Vp _ 2 \k are the piecewise polynomial functions of degree p — 2 that realize the 
bound (flTJl) with i = 0 on each K £Th- An easy adaptation of Lemma ITT! (and so also of Lemma 
4.1 in |T] or Lemma 3.1 in [7]) implies that, given p = max(l,p — 2), 


,min((p—2)+l,fc+l) 

(f,Vh) o.n - (fh,v h ) <C — K 


< C 


K&Th P‘ 

Lmin(p,/c+2) 


'k+1 


k+l,K p \ v h\l,K 


P 


k +2 


£ 

\KeT h 


k+l.K \ Vh ^ K - 


The final result follows by the definition of Th in (fl3l) and substituting p with p, up to a change 
of the constant C. □ 


By observing that, if the solution u of © is in H k+1 (UL) then f £ H k 1 (0), Lemma 14.41 
immediately gives also the following corollary. 

Corollary 4.5. Under assumptions (Dl) and (D2), let the solution u of (j2j be in H k+1 (Ll), 
k £ N. Then it holds 

_ h^ 

?h < C(fc, 7 , 7 )—r-||u|| fc+ i n, P = min(p, k). (27) 

pK 

Finally, we note that an analogous observation as in Remark [2] and Remark [3] holds also for 
Corollary 14.51 yielding 


Fh <C(k, 7 , j)-j-\u\k+i,n , l<k+l<p+l. (28) 

pK 

Remark 4. We stress the fact that using the same enhancing strategy introduced in [3j it is possible 
to obtain a more accurate load approximation. Nevertheless, the global order of convergence of 
the method does not change due to the presence of the other terms in the error estimates. 

4.2 hp estimate in the energy norm 

Finally, we are able to show the following convergence result. 

Theorem 4.6. Let k £ N, k > ^ and let the mesh assumptions (Dl) and (D2) hold. Let u and 
Uh be respectively the solution of problems ( 0 ) and m, with u £ H k+1 (Ll). Then, the following hp 
estimates hold: 


hn 

|u — Uh |i,n < Ci—rlMU+i.n, 

pK 

h k 

I u — Uh |i,n < C 2 — rMfc+i.a, 

pK 


p = min(p, k), 
ifk<p, 


(29) 

(30) 


where C\ and C 2 are two constants independent on h, p and u. 

Proof. It suffices to combine (fbll) . (fl3l) . (12T71) . (E21) and (177]l . □ 

Remark 5. Let the domain Q be convex. Following the argument shown in [10] (and, if p = 1,2 
suitably changing the definition of the discrete loading term (1111)1 and applying approximation 
results similar to those shown above, one can also easily derive L 2 estimates of the form: 

_ h 

||m — WfcJlo.fi < C(fc,7,7)— A t = min (p,k), (31) 

with the usual modification for the case k < p and where C is a constant independent on h and p. 





5 Exponential convergence for analytic solutions 

In this section, we derive an exponential convergence result for analytic solutions, under a further 
regularity assumption on the decomposition. We recall that we are given a polygonal decomposition 
Th and a triangular auxiliary subdecomposition Th (described in Section 14.1.21) . Given K polygon 
in T h ; we define Q = Q(K) as any of the smallest square containing K ; besides, given K triangle 
in Th, w e define Q = Q(K) the parallelogram given by Q = K UK*, where K* is the reflection 
of I< with respect to a midpoint of anyone of its edges. We point up that there are three possible 
Q(K)\ we fix arbitrarily one of them. Next, we define: 



We observe that dist(x, 12) < d(h), Vx £ 12 e xt, d(-) being a non-decreasing function in h. Therefore, 
\/h < h one has d(h) < d(h) and thus 12 ext is an uniformly bounded domain in terms of h, if h is 
bounded. We demand for the following regularity assumption on the mesh: 

(D3) there exists TV £ N independent on h such that there are at most TV overlapping squares 
in the collection {Q(K)} and TV parallelograms in the collection {Q(K)}, i.e. for all Q{K) 
in {Q(K)} and for all Q(K) in {Q(K)}, given Ik* '■= {Q(K) | Q(K) D Q(K') ^ 0} and 
I R , := {Q(K) | Q(K) D Q{K r ) ± 0}, it holds that card (I K ,), card(Ijj,)< N, VA £ T h and 
VKe%. 

We note that, given u £ H k+1 (Q ), k £ N, under assumption (D3), the following holds: 

H Mt+ 1 ,Q(K) < N\\u\\l +hnext , IMIfc+i,Q(jr) ^ JV 'll u llfc+i,n-,» 

K ^ T h KeT h 

with f2 ex t defined in (1321) . In order to obtain exponential convergence estimates for analytic func¬ 
tions, we must show bounds analogous to (OHl) and (IHD by expliciting the dependence of the 
constants Ci and C 2 on k, i.e. on the Sobolev regularity of the solution u. For this reason, we 
split this section in two parts. In Section [5711 we derive estimates of type (0H1) and (I5U1) with the 
dependence on k explicited; in Section [5.21 we derive an exponential convergence estimate. 

5.1 hp estimate using an overlapping square method 

In this Section, we use an overlapping square technique which allows us, under assumption (D3), 
to explicit the dependence on the Sobolev regularity in the estimate proven in Lemma 14.11 (con¬ 
sequently also on Lemmata 14.21 and 14.41) and on Lemma 14.31 Finally, we restate Theorem 14.61 on 
a proper extended domain, under assumption (D3). We note that the polynomial approximation 
which allows to have an estimate in p will be different from that discussed in Section [H such a 
polynomial, introduced by Babuska and Suri in [5] and [7], is a Fourier-type approximation. We 
decide to use here a different choice, by choosing an approximant of Legendre type whose properties 
are studied for instance in ESI■ The reason for this change is discussed in Remark [G] 

5.1.1 A first local estimate 

Here, we give an explicit representation of the constant C in (1151) in terms of k , k being the Sobolev 
regularity of the target function. We start by showing the counterpart of Lemma 14.11 As a minor 
note, we point out that the estimate of Lemma 15.11 does not require explicitely a shape regularity 
condition on the polygons, differently from Lemma 14.II 

Lemma 5.1. Let Q be the square [—1, l] 2 . Let K C Q be a polygon with barycenter ~x. R = 0. 
Moreover, assume that p > 2k, with k £ N. Then, there exists a family of projection operators 
{fi § },p=l,2,... with : H 2 (Q ) — > P p (Q) such that, for any u £ H k+1 (Q), it holds 

\ u ~^-q,p u \i,k ~ C2. e p "M k+i,Q (33) 

with C a constant independent on u, k and p. 
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Proof. Let Q = [—1, l] 2 . Let [Vi}f =1 be the set of vertices of Q. Let u £ H k+1 (Q). Let Q p (0) be 
the set of polynomials of maximum degree p in each variable over a domain 0 G R 2 . As an easy 
consequence of Lemma 4.67 in [2B], it is possible to show the existence of (p p £ Q P (Q) such that: 

PpW) = u(Vi), Vi = 1,... ,4 (34) 


and 


V’pIi,Q S 


< 2 


(p-k)l 


(p - k + 1 )! 


(p + k)\ p(p+ 1 ) (p+k — 1 )!J k + 1 ’Q 

Since p > k, it is easy to show that (1551) leads to the following simpler bound: 


(35) 


\u-(p p \ 1 Q < Cc k p k \u\ k+lt Q, with C = y/e. (36) 

In order to show this, we perform the computations only on the first term in the right hand side 
of (1551) since the treatment of the other one is analogous. Using Stirling fomula: 


( P~k)\ 
(p + k)\ 


(p + fc)(P +fe ) • e“(P+ fc ) • y/2n(p + k ) • e 9p + k 12 n + 1 12 n 


Then: 

< P~ 2k ■ e 2k ■ e 9p ~ k < Ce 2k p~ 2k , with C = e. (37) 

(p + k )! 

At this point, we observe that Q P (Q) Q P 2 p(Q)- This fact and (1551) immediately imply that there 
exists p p £ P P (Q) which interpolates u at the vertices of Q as in (1551) and which satisfies 

\u-p P \ h Q < C2 k e k p~ k \u\ k+1 Q , 

provided that p > 2k. We note that, owing to the fact that K C Q 1 it holds 

l M — ^Pp\l,K — \ U — Pp\l.Q — e P \ u \k+l,Q' 


In order to conclude, it suffices to define := tp p . □ 

The counterpart of Lemma l4~2l follows. 

Lemma 5.2. Let K £ Th- Let Q = Q(K) be the smallest square containing K and let u £ 
H k+1 (Q). Let p > 2k. Then, there exists a sequence of projection operators LIq p , p = 1,2,... 
with IIq p : H 2 (Q) — > P P (Q) such that for any k £ N 

h>* 

|«- n Q, p w|i,if < CM k -^-\\u\\ k+ i^ Q , p = min(p, k), 
where C and M are two constants independent on k, h, p and u. 

Proof. It suffices to apply Lemma |5.II and a classical scaling argument. The mapping F between 
Q and Q is the composition of a rototraslation and a dilatation in R 2 . The polygon K £ Q 
and the operator p u will be simply given by K = F(K ) and II q p u = (IIq (u o U -1 )) o F 
respectively. □ 

As done in Section H.l.ll we define u„ £ S k ’~ 1 (Th), Sh’~ 1 (Th) being introduced at the beginning 
of Section 14.11 as 

u„\k = {Uq ^k, with Q = Q(K), VK £ Th- 

Owing to assumption (D3) and Lemma U~?l we are able to give the following global estimate: 


h^ 

u-11^11,1,(1 < CA k —j:\\u\\k+i,n mt , p = mm(p,k), 


(38) 


where Ll ex t is defined in (1521) and C and A are two constants independent on h , p, k , 7 , 7 and u 
(A is independent also oniV). 
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5.1.2 A second local estimate 


In the present section, we give an explicit representation of the constant C in (1221) in terms of 
k. We point out that here the shape regularity assumption is needed; in fact, the usual scaling 
arguments used herein are based on affine mappings of shape regular triangles into the master 
triangle. 


Lemma 5.3. Under assumptions (D1), (D2) and (D3), provided that p > 2k, there exists ui £ Vh 
such that 

h k 

|« - w/|i,n < C ■ B k — \u\ k +i,n ext , (39) 

pK 

where C and B are two constants independent on k, p, h and u (B is independent also on N). 


Proof. The proof of this lemma is a combination of the arguments used in Lemma 15.11 and the 
construction of Lemma PPl Therefore, we only givejdrc sketch of the proof. We start by considering 
a triangle K in the subtriangular decomposition Th, we map it into the master triangle T (i.e. 
the triangle obtained halving the square [—1,1 ] 2 through its diagonal), we use a Legendre-type 
approximant in order to derive a estimate in p as in Lemma 15. 11 we go back to the triangle K. Let 
Q be the parallelogram Q = Q(K) (see assumption (D3)) and let {Vj }® =1 be the set of the vertices 
of K. Therefore, it is possible to show the existence of a £ P p (A) such that tpp(Vi) = u(Vi), 
Vi = 1,2,3 and such that 


- hl < CB k - 


\ u -^p\i,R S 


lfc+l,Q> 


(40) 


where C and B are two constants independent on p , h, k and u (B is also independent on N, 7 and 
7 ). We point out that this estimate holds for all the triangles in the triangular subdecomposition 
Th- We denote, with a little abuse of notation, by : fl —> R the global piecewise polynomial 
function whose restriction on each triangle K satisfies m- 

So far, we have obtained discontinuous piecewise polynomials. We set: 


E = E{K) 



\ 

U e 


^KeT h \KnK=ej 

/ 


U Q(K), 


e £ £ 


K ’ 


where we recall that is the set of the edges of K and Q(K) is defined in assumption (D3). We 
need to modify Lp k in order to get a continuous piecewise polynomial over Th without changing the 
approximation property cop . This can be done following the same approach as in El Theorem 
4.6, Lemma 4.7], i.e. by correcting ip k with suitable polynomial extensions of its edge jumps. It is 
easy to check that such step does not introduce constants depending on k. 

With another little abuse of notation, we have obtained a piecewise continuous 

polynomial of degree p over the subtriangular decomposition Th, such that an analogous of m 
holds for all K £%.'■ 

~ k h k 

|w- Pp\ ly K ^ c ( 7,7)5 prMfc+b-E- 

Using assumption (D3) and the arguments described in Lemma 14.31 it is easy to conclude the 
proof. □ 


The counterpart of Lemma 14.41 follows easily from Lemma 14.41 and Lemma 15.21 In particular 
the following holds. 

Lemma 5.4. Under assumptions (D1), (D2) and (D3), let f l ext be defined in (13211 . let the loading 
term f £ H k+1 (£l ex f). Then it holds 

1 h^ 

T h <CD~-^\\f\\ Ui ncxt , p = min(p, k + 2), (41) 

where C and D are two constants independent on k, h, p, 7 , 7 and u (D is also independent on 
N). 
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5.1.3 A global estimate result 

Combining bounds (1551) . (1551) . (TUI) and (1151) yields the following result. 

Theorem 5.5. Let k G N, k > -j. Let the mesh assumptions (Dl), (D2) and (D3) hold. Let u 
and Uh he respectively the solution of problems flj ) and Op. with u G H k (Tl). Let Ll ext he defined 
as in (1321) . Let u G H k+1 (LI ex f). Let 7, 7 and N be the constants introduced in assumptions (Dl), 
(D2) and (D3). Assume also p > 2k. Then, the following hp estimate holds: 

\u - u h |i,n < CA k — \u\ k +i,n ext , (42) 

pK 

where C and A are two constants independent on h, p, k and u (A is also independent on N). 

As done in Remark [5j we point out that if the domain Q is convex it is possible to derive easily, 
owing to the approximation properties of Legendre polynomials, L 2 estimates of the form: 

\u - u h \o,n < CA k —^\u\ k +i,n ext , (43) 

where C and A are two constants independent on h, p , k and u ( A is also independent on N ). 

Remark 6. We point out that in order to obtain the hp estimates of Theorem 14.61 and of Theorem 
15.51 we used two different approximant polynomials. Throughout Section [U we decided to follow 
the Babuska-Suri construction (see [B] and 0 ) which is based on a Fourier series expansion on a 
proper domain; this choice is essentially a matter of taste but has the merit of avoiding to use 
bi-polynomial functions. Nevertheless this construction obliges, also in the case of the overlapping 
square technique introduced at the beginning of Section [5j to use some extension operator (for 
instance the one described in m for Lipschitz domains). Thus, to give an explicit representation 
of the dependence of the involved constant on the Sobolev regularity k is a not trivial work. On 
the other hand, throughout Section [5j we made use of Legendre-type approximant (as done for 
instance in ESI)- In this case, owing to Legendre polynomials properties, we are able to obtain 
exponential estimates (see Section 15721) . since the dependence in the constant with respect to the 
Sobolev regularity k can be derived. 

5.2 Exponential convergence 

We have the following exponential convergence result for analytic solutions u over the extended 
domain fi ext (see (1551) 1. 

Theorem 5.6. Let the mesh assumptions (Dl), (D2) and (D3) hold. Letu anduh he respectively 
the solution of problems 0) and m, with u G A(flext), Aiytext) being the set of analytic function 
over the closure of Ll ext defined in m- Then, the following exponential convergence estimate 
holds: 

\\u-u h \\ hn <Ce- bp , (44) 

for some positive constants C and b independent on p. 

Proof. We recall (see for instance M) that an analytic function in the closure of a domain 0 G R 2 
is characterized by the following bound: 

||^“u|| 00i e<CAl“la!, a = ( ai ,a 2 )GN 0 2 , (45) 

where a\ = aqla^! and where C and A are constants independent on the multi-index a; neverthe¬ 
less, C and A depends on u and on 0. Recalling (1551) . we have 

~ h k 

\u-u h \ h n < (7(7,7, A)A( 7 , 7 ) fc - F |u| fc+ i,n ext , 
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if p > 2k. Using standard results from space interpolation theory [13][28|, from the above bound 
one can easily derive 

\u-u h \i,n < C(7,7,7V)A(7,7) s ^ 7 |u| s+ i ) n ext (46) 

pS 

for all 5 G R with 2 < 2 s < p. The combination of Si) and © yields 


— Uh\i,n < C 41 S+I (s + 1)!. 


By means of Stirling formula, we obtain: 

(haaV f s + iy +1 


^ / hAA\ fs+ lV +i r—, 

u - w+.n < C I —— J (— — J V / 27r(s + l) 2 , 


easily yielding 


I , . r hAA' 

I u - u h < 6 - s 

V e P ; 


By denoting 5 = we can write: 


I u- u h |i,n < C ( -S ) s*. 

,P 


Since this last inequality holds true for all s such that 2 < 2 s < p, we may choose s 
Hence: 


|u - u h |i,a < C 


2(5+1) 


2 ( 5 + 1 ) 


3 

P 2 


Ce~ bp pi , with b = 


10 S(2(£t)) 
2(5 + 1) 


— R 
2 ( 5 + 1 )’ 


(47) 


The multiplier pi can be absorbed by e bp by making b a little bit smaller and increasing C; 
therefore, © immediately yields 

\u-u h \ h n<Ce~ bp , (48) 

for some constants C and b independent on p. The result follows by the Poincare inequality. □ 


6 Numerical results 

In this section, we present numerical results experimentally validating the error estimates ®, 
©, © and (14811 . We consider four types of meshes (see Figure [1]) on the domain Q = [0,1] 2 , 
namely an unstructured triangular mesh, a regular square mesh, a regular hexagonal mesh and 
a Voronoi-Lloyd mesh (see [H]). The basis {<? a } of the space P p _ 2 (iF) introduced in Section [3] 



Figure 1: From left to right: unstructured triangular mesh, regular square mesh, regular hexagonal mesh, Voronoi- 
Lloyd mesh 


\/K G Thi is taken to be the same as that introduced for instance in [5] or m- Different choices 
are investigated in the Appendix. Moreover, we fix a possible choice for the stabilizing term S K 
introduced in m (see for instance 0 ) as: 

dim (Vh(K)) 

S K (u h ,v h ) = ^2 Xr(u h )xr(v h ), Vu h , v h e V h (K), K e +, (49) 

r=l 
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where \r , Vr = 1,...,dim(V/j,(if)) is the operator which associates to each function in the local 
space Vh{K) its rth local degree of freedom. 

Remark 7. The stabilizing bilinear form does not guarantee the stability property (fTOl) uniform 
in p. Nevertheless, the numerical results seem to be robust with respect to this choice. A theoretical 
study of the stabilization will be the object of further investigation. 

In order to estimate the error introduced by the MATLAB algebraic sparse solver, we have 
solved a problem whose exact solution is the polynomial u(x, y) = x 2 + y 2 . Since the VEM passes 
the patch test, in this case, for k > 2, the approximate solution Uh coincides with u and the error 
that we measure is only due to the algebraic solver (it should be zero in exact arithmetic). Hence, 
together with the standard error in the ff-norm and in the L 2 -norm with respect to the solutions 
(1501) and m , in our convergence figures we also plot this algebraic error. When the error curve 
comes close to the algebraic error curve, the convergence error and the error introduced by the 
MATLAB algebraic solver are of the same order and the expected theoretical behaviour does not 
hold anymore. 

6.1 Convergence in p for an analytic function 

We consider problem © with loading term f(x,y) = 2n 2 sm(irx)sm(ny). The exact solution is 
given by: 

u(x, y) = sin(7ra;) sin(7ry). (50) 

In this tests, the mesh is kept fixed (see Figure © and the polynomial degree is raised. In Figures 
[2] and [3] we report the errors among the discrete and exact solutions. Since we are dealing with a 
virtual element solution Uh (that is unknown inside elements) we cannot direcly compute the error 
\\u — u/j||s, 12, s = 0,1. Therefore, as is standard in VEM, we plot instead ||u — nJu^Ho.n and 
|u — that are good representatives of the above errors (see © for the definition of the 

operator nj and (fl^l) for the definition of the H 1 broken Sobolev seminorm). 




Figure 2: u(x,y) = sin(Trr) sin(Try): unstructured triangle mesh (left); regular square mesh (right) 


In accordance with Theorem 15.61 the exponential convergence is evident from the decreasing 
slope in the error graphs. Moreover, from Figures [2] and [3] we can observe that the lower line is a 
good marker for the indication of the machine algebra error. 

6.2 Convergence in p for a function with finite Sobolev regularity 

Secondly, we present a similar behaviour test for the case of a problem with solution 

u(r, 6) = r 2 ' 5 sin(2.50), (51) 

where (r, 6) are the polar coordinates with respect to the origin. Since the function is harmonic, 
the loading term f = 0 and the Dirichlet boundary conditions are set in accordance with u\on- 
We note that u £ H 3 - 5 ~ E (Q), Ver > 0. In Figures 0] and [5j the segmented line represents a line of 
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Figure 3: u(x,y) = sin(7nr) sin(7ry); regular hexagonal mesh (left); 4: Voronoi-Lloyd mesh (right) 


slope 5 = 2- 2.5. Owing to Theorem 14.61 we should have an estimate in p of the type p ~ a , a = 2.5. 
Anyhow, for this type of corner singularity, one could extend the technical result of [B] and [7] 
obtaining error estimate in p of the type p~ 2a also in our VEM framework. Thus, we expect a 
slope for the H 1 error of the type p ~ 5 , which is represented with the dashed line in Figures 0] and 
[5| Figures 0] and 01 are in agreement with such observation. 




Figure 4: u = r 2,5 sin(2.50); unstructured triangle mesh (left); 2: regular square mesh (right) 


6.3 Convergence in h 

In this subsection, we show the convergence rate when the polynomial degree is kept fixed and 
the meshsize goes to zero. We consider a sequence of hexagonal and Voronoi-Lloyd meshes and 
we study the same harmonic test case as in Section 16.II In particular, we examine the case p = 3 
and p = 5. We observe that the slope of the errors are in accordance with Theorem 14.61 and with 
estimate (ED- The same considerations about Figure 0] are still valid for Figure 0] We only point 
out that the strange L 2 error behaviour for the final step is due to the machine precision error. 


7 Appendix 

In this appendix, we study the behaviour of the condition number of the global stiffness matrix 
of problem (0]) and we explore some alternatives in the choice of the local VEM basis. We note 
that, as it happens in Finite Elements (see 0] and 0), the main responsible for the growth of the 
condition number (when p increases) are the internal “bubble” basis functions. In Section 17.II we 
numerically investigate the behaviour of the condition number by changing the polynomial basis 
{q a } of P p _ 2 (AT) (for all K £ Th), introduced in Section 0] for the definition of the local internal 
degrees of freedom. In Section [DD we discuss an orthogonalization technique that strongly reduces 
the condition number, but is unstable with respect to machine precision. 
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Figure 5: u = r 2-5 sin(2.50); regular hexagonal mesh (left); 4: Voronoi-Lloyd mesh (right) 


io~ 2 


io~ 3 



io° io' io 2 

1/diam 



10°‘ 3 10°' 5 10°' 7 10° 9 
1/diam 


Figure 6: regular hexagonal mesh (left); Voronoi-Lloyd mesh (right); p=3 


7.1 Three explicit bases 


In this subsection, we consider three explicit basis: 

• {gi}, the same basis introduced for instance in [5] and mi. that is to say: 

q ^ = {~h~ S ) ’ Va € Ng, |a| < p- 2, (52) 

where xk and hx are respectively the barycenter and the diameter of the polygon K. 

• {Qa }, which is defined by 


q 


2 

CX. 


u 


|0 ,K 


Va € Nq, |a| < p — 2. 


(53) 


{(Jab which is a Legendre-type basis. In order to define it, we recall that Nk is the number 
of vertices of K and { V % }is the set of vertices of K; moreover, we set 

- ~ . I Xy ,max Xy m i n MV. max 

Xk = (X K ,yK) = 


— |^V,max + ^'V,min|; ^x — | //V.rnax T UV , 


^N K 


Nk . 


,N K , 


where zy jmax = max^ij, xy min = irim^ x t , yy m ax = max:^^, y v , min = min X=\Vi- 
Besides, let L s (-) be the Legendre polynomial of degree s on the segment [—1,1] (see e.g. [33] 
for the properties of Legendre polynomials). Then, we are able to define the basis: 


Qa — 


x — Xx 


l K 


L n 


,y-yK 

' h v 

n K 


Mol € Nq, |a| < p — 2. 


(54) 
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Figure 7: regular hexagonal mesh (left); VoronoiLloyd mesh (right); p=5 


We observe that the third choice should be, at least at a first glance, better than the other two, 
thanks to orthogonality properties of Legendre polynomials. We will see that instead this is not the 
case in general. Indeed, the orthogonality properties of the Legendre basis are quickly lost when 
the considered domain is not rectangular. In our tests, we consider the same meshes already used 
in Section [HI see Figure [TJ In Figures [8] and GO we compare the behaviour of the condition number, 
given by the Matlab command cond, for the three choice of the bases mentioned above and the four 
meshes. We stress the fact that the Legendre-type basis performs better in the case of the square 




Figure 8: 1: unstructured triangle mesh; 2: regular square mesh 


mesh; this is believable thanks to the orthogonality properties of the Legendre polynomials. On 
the contrary, one can see that more general meshes, such as the hexagonal, unstructured triangular 
and Voronoi-Lloyd ones, the best result is obtained with the L 2 -scaled basis. Finally, we present 
in Figures [TOJ [TT] fl2l and flTTI some “comparison tests” in which we report the error \u — IlJ Uh |i,n 
(see ©), by using the four different meshes in Figure [T] and the three different bases. We consider 
the same two test cases of Section 0 From Figures [TO] and [TT] one can see that the Legendre- 
type basis is a good choice for the square case, while on triangles it is very unstable; besides, for 
general meshes, it seems that the slope of the error with the other two bases is almost the same 
and performs better than the Legendre-type basis. The same considerations for Figures [TU] and 
[Tl]hold for Figures ITT and ITT! 

We point out that in our numerical tests we have used a direct solver in order to work out the 
global system arising from the discrete problem ([3]). A consequence of this fact is that the condition 
number of the global matrix does not affect the resolution of the linear system as it would do if we 
used an iterative solver. This explains why the behavior of the errors with the choice of the classical 
basis { q and the scaled basis {g^} is almost the same, notwithstanding the large difference in 
the condition number of the global matrix as shown in Figure [8] and in Figure 0 
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Figure 9: 3: regular hexagonal mesh; 4: VoronoiLloyd mesh 




Figure 10: u(x,y) = sin(7nc) sin(7r?/); 1: unstructured triangle mesh; 2: regular square mesh 


7.2 A “Virtual” Gram Schmidt process 

From the previous subsection, it is clear that by suitably changing the basis of the space one 
obtains better condition numbers for the global stiffness matrix. Despite this, from Figures [5] and 
M we note that the condition numbers are still large and, although in our codes we use a direct 
solver, it would be preferable to reduce such numbers. Therefore, in this subsection, we consider a 
technique which considerably reduces the condition number of the global stiffness matrix, but at 
the price of an unstable propagation of the machine error precision (as better discussed later). In 
particular, we begin working locally and we fix an element K G Tk] we consider the local stiffness 
matrix . We recall that the th entry of A£ has the following form (see [9]): 

[Ah]i,j ■= [A ]ij = afiiipuVj) = a K {HJ (pj) + S K {<pi - Ft - Ft Jpj). 

Moreover, we recall that 


g._ n B d + I d ,B d , + I d 

where Id and Bd are respectively the internal and the boundary degrees of freedom. As already 
mentioned, the internal “bubble” basis functions are the main responsible for the growth of the 
condition number when p increases. 

In Section [7.2.11 we introduce a technique which allows us to orthonormalize the internal elements 
of the local basis with respect to the discrete bilinear form ■ In Section 17.2.21 we present 
numerical experiments in which we compare the condition numbers and H 1 type errors of both 
standard VEM (i.e. those using the basis {g^}) and this new choice. 
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Figure 11: u(x,y) = sin(Tr.x) sin (Try); 3: regular hexagonal mesh; 4: VoronoiLloyd mesh 




Figure 12: u(r,9) = r 2 5 sin(2.50); l:unstructured triangle mesh; 2: regular square mesh 


7.2.1 A “virtual” Gram Schmidt method 

Let Gd , Bd and Id be the number of total, boundary-vertex and internal degrees of freedom 
respectively over K of the space 14(A). Let 

= Wf }f=r U {YiYiti (55) 

be the basis of 14(A) introduced in [IT] , where i represents the “boundary” basis elements 

and {iPiYiL\ represents the “internal” basis elements. We want here to introduce a new basis which 
shares the same “boundary” elements of (l55l) and which has orthonormalized “internal” elements 
with respect to the discrete local bilinear form ■ The new basis reads 

{&}£?i = Wj}f= 1 u {<Pi Yit 1> (56) 

with the <pj are to be defined. Since we want ortho normality with respect to , we simply apply 
the Gram Schmidt method to the internal elements of (l55l) . We observe that each Ipf can be split 
as: 


<Pl — 

k=1 

where the Aare the Gram Schmidt coefficients. 

Before giving an explicit representation of the A k,i, we want to observe the following fact. Let A// 
be the matrix having (*, j)-th entry given by (Jp\, <Pj). Then, thanks to the choice of using the 

Gram Schmidt method, we have that A//, is a diagonal matrix. We can explicitly represent such 
matrix by means of: 

A a = A • An • A 1 , 
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Figure 13: u(r, 0) = r 2 5 sin(2.50); 3: regular hexagonal mesh; 4: VoronoiLloyd mesh 


where A is the left lower triangular matrix having in the l -th row the coefficients k == 1,..., l 

of <?{■ 

Out task, now, is to find the matrix A. We stress the fact that in order to compute such matrix we 
can not use explicitly the “internal” basis elements since they are virtual; notwithstanding, we can 
only use the scalar products We start by orthogonalizing the basis. We know that: 


i-1 


H>\ = v>{, 


<Pl 






(<Pi 




k—1 




\/i 


2 


By induction, one can show that the coefficients of the function reads: 


Az,z = 1, 


A; 


i -1 

= - H 

j—k 


“hWi'Pj) 
°h ($>$)’ 


Vk = 1,..., Z — 1. 


(57) 


We observe that the coefficients Xj t k in m are already known at the 1-th step. On the other 
hand, we need a routine which permits to compute ajf (ipj. <p T j ) and , £>f), without knowing 

explicitly the basis and only knowing the scalar products One has: 


a h 


j j 

= a h(J 2 x j,hVh) = A Cb 0 • A ii' A T (:,j) 

i— 1 h= 1 


(58) 


and 


-.K 


{<Pl,<Pj) = a h{v< 


l 7 


= An(2,:) • a T (:, j). 


(59) 


Thus, we have the tools for finding out the coefficients of A. We remark that we have obtained 
a new local virtual basis and that also the internal moments in the definition of the degrees of 
freedom are changed. In particular, the basis of scaled monomials {m a } has changed to a new 
basis {to q }, but for the computation of the new local stiffness matrix, we do not need to know 
explicitly such a new basis. In fact, let I be the identity matrix having size equal to Bd. Then, we 
define: 


A = 


I 

0 T 


0 

A 


and the new local stiffness matrix is obtained by: 


A := 


A ss 

A ib 


A B / 

A H 


= A 


A bb 

A//3 


a BI 

A n 


•A t 


A BB A-bi 

A/b D 


It remains the normalization of the internal basis elements. We observe that the entries of the 
diagonal matrix Ajj are [A//],;,,; = (ipi, ipi). Therefore, we just need to multiply on the left and 
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on the right the matrix A by 


~ _ r 1 o 

k |A//|-iJ ’ 

where |A//| _ 5 is the diagonal matrix such that [A// I %]i,i = \J\A n \i,il obtaining 

~ _ r 1 

[o T A|A/j| _ 5_ ■ 

Finally, the definitive new local stiffness matrix reads: 


A = 

A bb 

A B1 

= A 

A bb 

A-bi 

>n 

II 

Abb 

A BI 


.A/-/j 

A n_ 


A-ib 

A h 


.A-ib 



( 60 ) 


( 61 ) 


It is clear that the global virtual basis is obtained with a classical conforming coupling of the 
boundary degrees of freedom. Finally, we stress the fact that in order to compute the projection 
nj, one may use the procedure described above; by means of straightforward computations, it is 
easy to recover analogous results to those presented in m- 


7.2.2 Numerical results using the “virtual” Gram Schmidt method 

In the present section, we present some numerical experiments on the behaviour of the condition 
number. We consider the four different meshes in Figure [l] and we compare the condition number 
of the L 2 scaled basis introduced in Section 17.11 and the new Gram-Schmidt basis, described in 
Section [7. 2 . II We remind that, from the previous numerical experiments, the L 2 scaled basis seems 
to be the most well conditioned among the choices of Section 17. 2. II From the results in Figures HUH 




Figure 14: 1:unstructured triangle mesh; 2: regular square mesh 


and it follows that the Gram Schmidt basis performs much better, at least for what concerns 
the condition number. 

In Figures [THl and flTl we compare the behaviour of the error \u— IlJ Uh |i.n using the two bases above 
on the usual test case u{x , y) = sin(7rx) sin(7ry). We observe that, although the method described 
in this subsection improves the condition number of the global stiffness matrix, it is numerically 
unstable. Therefore, in practice, the proposed Gram-Schmidt method may be preferable to the 
simple basis choice in (15611 only for mid-low values of p. 
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